Custom Figures for TS analysis
Javier Rodriguez Casariego
November 13, 2024
DE conditional upset plot
# Pull and organize DE data for all timepoints.
# For all timepoints (Ctrol, hyp_2h ...), contrast is always Lab_cross vs Wild_cross (upregulation means higher in lab_cross and viceversa)
knitr::opts_chunk$set(root.dir = '/Volumes/JARC_2T/NSF-CREST_Postdoc/LongTerm_Hypoxia_exp/AplCal_hypoxia')
library(ggupset)
library(ComplexHeatmap)## Loading required package: grid## ========================================
## ComplexHeatmap version 2.21.1
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors# load gene expression results for each timepoint for Lab_cross vs Wild_Cross
load("gene_expression/output/TS_output/LCvsWC_Abdominal/timeSeries_obj_LCvsWC_abdominal.Rdata")## load("gene_expression/output/TS_output/LCvsWC_Pleural/timeSeries_obj_LCvsWC_pleural.Rdata")
names(A_TS_object@DE_results$conditional)## [1] "Lab_cross_vs_Wild_cross_TP_0" "Lab_cross_vs_Wild_cross_TP_1"
## [3] "Lab_cross_vs_Wild_cross_TP_2" "Lab_cross_vs_Wild_cross_TP_4"
## [5] "Lab_cross_vs_Wild_cross_TP_5" "Lab_cross_vs_Wild_cross_TP_6"
## [7] "Lab_cross_vs_Wild_cross_TP_3"names(P_TS_object@DE_results$conditional)## [1] "Lab_cross_vs_Wild_cross_TP_0" "Lab_cross_vs_Wild_cross_TP_1"
## [3] "Lab_cross_vs_Wild_cross_TP_2" "Lab_cross_vs_Wild_cross_TP_4"
## [5] "Lab_cross_vs_Wild_cross_TP_5" "Lab_cross_vs_Wild_cross_TP_6"
## [7] "Lab_cross_vs_Wild_cross_TP_3"Abdominal_Control_T0_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_0"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Control_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
Abdominal_Hyp_2h_T1_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_1"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_2h_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
Abdominal_Hyp_6h_T2_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_2"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_6h_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
Abdominal_Reox_T3_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_3"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Reoxygenation_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
Abdominal_Hyp_5d_T4_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_4"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_5d_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
Abdominal_Hyp_6d_T5_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_5"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_6d_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
Abdominal_Recovery_T6_DE <- A_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_6"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Recovery_Abdominal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
# Pleural/Pedal
PleuralPedal_Control_T0_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_0"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Control_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
PleuralPedal_Hyp_2h_T1_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_1"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_2h_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
PleuralPedal_Hyp_6h_T2_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_2"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_6h_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
PleuralPedal_Reox_T3_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_3"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Reoxygenation_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
PleuralPedal_Hyp_5d_T4_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_4"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_5d_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
PleuralPedal_Hyp_6d_T5_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_5"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Hypoxia_6d_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)
PleuralPedal_Recovery_T6_DE <- P_TS_object@DE_results[["conditional"]][["Lab_cross_vs_Wild_cross_TP_6"]]$DE_raw_data[,c('gene_id','log2FoldChange','padj')] %>%
filter(.,padj < 0.05) %>%
mutate(dir = if_else(log2FoldChange > 0, "UP", "DOWN")) %>%
mutate(contrast = paste0("Recovery_PleuralPedal","_",dir)) %>%
dplyr::select(gene_id, contrast, dir)#Abdominal
# Create list with DE genes per contrast
DE_lt_A <- list(Normoxia = as.character(Abdominal_Control_T0_DE$gene_id),
Hypoxia_2h = as.character(Abdominal_Hyp_2h_T1_DE$gene_id),
Hypoxia_6h = as.character(Abdominal_Hyp_6h_T2_DE$gene_id),
Hypoxia_5d = as.character(Abdominal_Hyp_5d_T4_DE$gene_id),
Hypoxia_6d = as.character(Abdominal_Hyp_6d_T5_DE$gene_id),
Reoxygenation = as.character(Abdominal_Reox_T3_DE$gene_id),
Recovery = as.character(Abdominal_Recovery_T6_DE$gene_id))
DE_lt_A_down <- list(Normoxia = as.character(Abdominal_Control_T0_DE$gene_id[Abdominal_Control_T0_DE$dir == "DOWN"]),
Hypoxia_2h = as.character(Abdominal_Hyp_2h_T1_DE$gene_id[Abdominal_Hyp_2h_T1_DE$dir == "DOWN"]),
Hypoxia_6h = as.character(Abdominal_Hyp_6h_T2_DE$gene_id[Abdominal_Hyp_6h_T2_DE$dir == "DOWN"]),
Hypoxia_5d = as.character(Abdominal_Hyp_5d_T4_DE$gene_id[Abdominal_Hyp_5d_T4_DE$dir == "DOWN"]),
Hypoxia_6d = as.character(Abdominal_Hyp_6d_T5_DE$gene_id)[Abdominal_Hyp_6d_T5_DE$dir == "DOWN"],
Reoxygenation = as.character(Abdominal_Reox_T3_DE$gene_id[Abdominal_Reox_T3_DE$dir == "DOWN"]),
Recovery = as.character(Abdominal_Recovery_T6_DE$gene_id[Abdominal_Recovery_T6_DE$dir == "DOWN"]))
DE_lt_A_up <- list(Normoxia = as.character(Abdominal_Control_T0_DE$gene_id[Abdominal_Control_T0_DE$dir == "UP"]),
Hypoxia_2h = as.character(Abdominal_Hyp_2h_T1_DE$gene_id[Abdominal_Hyp_2h_T1_DE$dir == "UP"]),
Hypoxia_6h = as.character(Abdominal_Hyp_6h_T2_DE$gene_id[Abdominal_Hyp_6h_T2_DE$dir == "UP"]),
Hypoxia_5d = as.character(Abdominal_Hyp_5d_T4_DE$gene_id[Abdominal_Hyp_5d_T4_DE$dir == "UP"]),
Hypoxia_6d = as.character(Abdominal_Hyp_6d_T5_DE$gene_id)[Abdominal_Hyp_6d_T5_DE$dir == "UP"],
Reoxygenation = as.character(Abdominal_Reox_T3_DE$gene_id[Abdominal_Reox_T3_DE$dir == "UP"]),
Recovery = as.character(Abdominal_Recovery_T6_DE$gene_id[Abdominal_Recovery_T6_DE$dir == "UP"]))
#calculate the intersection between the differentially expressed gene sets
set.seed(12345)
DE_m <- make_comb_mat(DE_lt_A, mode = "distinct")
DE_m_down <- make_comb_mat(DE_lt_A_down, mode = "distinct")
DE_m_up <- make_comb_mat(DE_lt_A_up, mode = "distinct")
# exclude self-intersects (total # of diff. genes will be displayed separately)
DE_m <- DE_m[comb_degree(DE_m) > 1]
DE_m_down <- DE_m_down[comb_degree(DE_m_down) > 1]
DE_m_up <- DE_m_up[comb_degree(DE_m_up) > 1]
DE_m <- DE_m[comb_size(DE_m) > 20]
DE_m_down <- DE_m_down[comb_name(DE_m)]
DE_m_up <- DE_m_up[comb_name(DE_m)]
m_list <- list(all_DEGs = DE_m,
downregulated = DE_m_down,
upregulated = DE_m_up)
sapply(m_list, comb_size)## all_DEGs downregulated upregulated
## 1111111 162 116 46
## 1111011 59 19 40
## 0111111 297 271 26
## 1110011 33 2 31
## 0111011 152 107 45
## 0011111 22 17 5
## 0111010 21 14 7
## 0111001 60 41 19
## 0110011 58 19 39
## 0011011 52 32 20
## 0110010 41 18 22
## 0110001 97 15 82
## 0101001 32 14 18
## 0100011 26 12 14
## 0011001 43 32 9
## 0010011 45 25 20
## 0001011 48 16 32
## 0110000 127 51 77
## 0101000 31 13 17
## 0100010 33 16 15
## 0100001 98 22 67
## 0011000 43 21 18
## 0010010 36 22 14
## 0010001 89 48 41
## 0001010 31 14 16
## 0001001 92 44 49
## 0000011 57 31 25m_list = normalize_comb_mat(m_list)
sapply(m_list, comb_size)## all_DEGs downregulated upregulated
## 1111111 162 116 46
## 1111011 297 271 26
## 1101111 59 19 40
## 1101011 152 107 45
## 1001111 33 2 31
## 0111011 22 17 5
## 1101010 60 41 19
## 1101001 21 14 7
## 1001011 58 19 39
## 0101011 52 32 20
## 1100010 32 14 18
## 1001010 97 15 82
## 1001001 41 18 22
## 1000011 26 12 14
## 0101010 43 32 9
## 0100011 48 16 32
## 0001011 45 25 20
## 1100000 31 13 17
## 1001000 127 51 77
## 1000010 98 22 67
## 1000001 33 16 15
## 0101000 43 21 18
## 0100010 92 44 49
## 0100001 31 14 16
## 0001010 89 48 41
## 0001001 36 22 14
## 0000011 57 31 25max_set_size = max(sapply(m_list, set_size))
max_comb_size = max(sapply(m_list, comb_size))
top_ha = HeatmapAnnotation(
"up" = anno_barplot(comb_size(m_list[[3]]),
gp = gpar(fill = "#B31B21"), height = unit(2.5, "cm"), ylim = c(0,350)),
"down" = anno_barplot(comb_size(m_list[[2]]),
gp = gpar(fill = "#1465AC"), height = unit(2.5, "cm"), ylim = c(0,350)),
"all" = anno_barplot(comb_size(m_list[[1]]),
gp = gpar(fill = "black"), height = unit(2.5, "cm"),ylim = c(0,350)),
gap = unit(1, "mm"), annotation_name_side = "left", annotation_name_rot = 90)
# the same for using m2 or m3
ht_A <- UpSet(m_list[[1]], top_annotation = top_ha,
right_annotation = upset_right_annotation(m_list[[1]], add_numbers = TRUE, show_annotation_name = F, axis = F),
set_order = c("Normoxia","Hypoxia_2h","Hypoxia_6h","Hypoxia_5d",
"Hypoxia_6d","Reoxygenation","Recovery"));ht_A
# Plot upset for all shared DE and directionality
pdf(file="gene_expression/figures/Fig2A_UpsetPlotDEGs_Abdominal.pdf", width = 8, height = 5)
draw(ht_A)
dev.off()## quartz_off_screen
## 2#Pleural_Pedal
DE_lt_P <- list(Normoxia = as.character(PleuralPedal_Control_T0_DE$gene_id),
Hypoxia_2h = as.character(PleuralPedal_Hyp_2h_T1_DE$gene_id),
Hypoxia_6h = as.character(PleuralPedal_Hyp_6h_T2_DE$gene_id),
Hypoxia_5d = as.character(PleuralPedal_Hyp_5d_T4_DE$gene_id),
Hypoxia_6d = as.character(PleuralPedal_Hyp_6d_T5_DE$gene_id),
Reoxygenation = as.character(PleuralPedal_Reox_T3_DE$gene_id),
Recovery = as.character(PleuralPedal_Recovery_T6_DE$gene_id))
DE_lt_P_down <- list(Normoxia = as.character(PleuralPedal_Control_T0_DE$gene_id[PleuralPedal_Control_T0_DE$dir == "DOWN"]),
Hypoxia_2h = as.character(PleuralPedal_Hyp_2h_T1_DE$gene_id[PleuralPedal_Hyp_2h_T1_DE$dir == "DOWN"]),
Hypoxia_6h = as.character(PleuralPedal_Hyp_6h_T2_DE$gene_id[PleuralPedal_Hyp_6h_T2_DE$dir == "DOWN"]),
Hypoxia_5d = as.character(PleuralPedal_Hyp_5d_T4_DE$gene_id[PleuralPedal_Hyp_5d_T4_DE$dir == "DOWN"]),
Hypoxia_6d = as.character(PleuralPedal_Hyp_6d_T5_DE$gene_id)[PleuralPedal_Hyp_6d_T5_DE$dir == "DOWN"],
Reoxygenation = as.character(PleuralPedal_Reox_T3_DE$gene_id[PleuralPedal_Reox_T3_DE$dir == "DOWN"]),
Recovery = as.character(PleuralPedal_Recovery_T6_DE$gene_id[PleuralPedal_Recovery_T6_DE$dir == "DOWN"]))
DE_lt_P_up <- list(Normoxia = as.character(PleuralPedal_Control_T0_DE$gene_id[PleuralPedal_Control_T0_DE$dir == "UP"]),
Hypoxia_2h = as.character(PleuralPedal_Hyp_2h_T1_DE$gene_id[PleuralPedal_Hyp_2h_T1_DE$dir == "UP"]),
Hypoxia_6h = as.character(PleuralPedal_Hyp_6h_T2_DE$gene_id[PleuralPedal_Hyp_6h_T2_DE$dir == "UP"]),
Hypoxia_5d = as.character(PleuralPedal_Hyp_5d_T4_DE$gene_id[PleuralPedal_Hyp_5d_T4_DE$dir == "UP"]),
Hypoxia_6d = as.character(PleuralPedal_Hyp_6d_T5_DE$gene_id)[PleuralPedal_Hyp_6d_T5_DE$dir == "UP"],
Reoxygenation = as.character(PleuralPedal_Reox_T3_DE$gene_id[PleuralPedal_Reox_T3_DE$dir == "UP"]),
Recovery = as.character(PleuralPedal_Recovery_T6_DE$gene_id[PleuralPedal_Recovery_T6_DE$dir == "UP"]))
#calculate the intersection between the differentially expressed gene sets
set.seed(12345)
DE_m <- make_comb_mat(DE_lt_P, mode = "distinct")
DE_m_down <- make_comb_mat(DE_lt_P_down, mode = "distinct")
DE_m_up <- make_comb_mat(DE_lt_P_up, mode = "distinct")
# exclude self-intersects (total # of diff. genes will be displayed separately)
DE_m <- DE_m[comb_degree(DE_m) > 1]
DE_m_down <- DE_m_down[comb_degree(DE_m_down) > 1]
DE_m_up <- DE_m_up[comb_degree(DE_m_up) > 1]
DE_m <- DE_m[comb_size(DE_m) > 20]
DE_m_down <- DE_m_down[comb_name(DE_m)]
DE_m_up <- DE_m_up[comb_name(DE_m)]
m_list <- list(all_DEGs = DE_m,
downregulated = DE_m_down,
upregulated = DE_m_up)
sapply(m_list, comb_size)## all_DEGs downregulated upregulated
## 1111111 311 248 63
## 1111101 21 17 4
## 1111011 46 11 35
## 0111111 83 72 11
## 1111001 22 4 18
## 0111011 74 32 42
## 0111010 29 15 14
## 0111001 48 13 35
## 0110011 55 25 30
## 0011011 22 8 14
## 0110010 63 18 45
## 0110001 70 17 51
## 0100011 23 12 9
## 0011001 42 19 23
## 0010011 43 19 23
## 0110000 105 38 69
## 0100010 21 13 10
## 0100001 59 28 26
## 0011000 34 15 19
## 0010010 84 24 61
## 0010001 67 38 27
## 0001001 57 24 33
## 0000011 36 20 16m_list = normalize_comb_mat(m_list)
sapply(m_list, comb_size)## all_DEGs downregulated upregulated
## 1111111 311 248 63
## 1111110 21 17 4
## 1111011 83 72 11
## 1101111 46 11 35
## 1101110 22 4 18
## 1101011 74 32 42
## 1101010 48 13 35
## 1101001 29 15 14
## 1001011 55 25 30
## 0101011 22 8 14
## 1001010 70 17 51
## 1001001 63 18 45
## 1000011 23 12 9
## 0101010 42 19 23
## 0001011 43 19 23
## 1001000 105 38 69
## 1000010 59 28 26
## 1000001 21 13 10
## 0101000 34 15 19
## 0100010 57 24 33
## 0001010 67 38 27
## 0001001 84 24 61
## 0000011 36 20 16max_set_size = max(sapply(m_list, set_size))
max_comb_size = max(sapply(m_list, comb_size))
top_ha = HeatmapAnnotation(
"up" = anno_barplot(comb_size(m_list[[3]]),
gp = gpar(fill = "#B31B21"), height = unit(2.5, "cm"), ylim = c(0,350)),
"down" = anno_barplot(comb_size(m_list[[2]]),
gp = gpar(fill = "#1465AC"), height = unit(2.5, "cm"), ylim = c(0,350)),
"all" = anno_barplot(comb_size(m_list[[1]]),
gp = gpar(fill = "black"), height = unit(2.5, "cm"),ylim = c(0,350)),
gap = unit(1, "mm"), annotation_name_side = "left", annotation_name_rot = 90)
# the same for using m2 or m3
ht_P <- UpSet(m_list[[1]], top_annotation = top_ha,
right_annotation = upset_right_annotation(m_list[[1]], add_numbers = TRUE, show_annotation_name = F, axis = F),
set_order = c("Normoxia","Hypoxia_2h","Hypoxia_6h","Hypoxia_5d",
"Hypoxia_6d","Reoxygenation","Recovery"));ht_P
pdf(file="gene_expression/figures/Fig2B_UpsetPlotDEGs_Pleural.pdf", width = 8, height = 5)
draw(ht_P)
dev.off()## quartz_off_screen
## 2Functional Enrichment Combined Figures
GO:BP Abdominal
library(clusterProfiler)
library(ggplot2)
library(dplyr)
library(stringr)
library(forcats)
library(GO.db)
# GO:BP Abdominal
GO_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseGO_results.tab")
GO_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseGO_results.tab")
GO_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseGO_results.tab")
GO_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseGO_results.tab")
GO_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseGO_results.tab")
GO_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseGO_results.tab")
GO_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseGO_results.tab")
gene_count.GO_T0 <- GO_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T1 <- GO_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T2 <- GO_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T3 <- GO_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T4 <- GO_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T5 <- GO_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T6 <- GO_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(GO_res_T0, gene_count.GO_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(GO_res_T1, gene_count.GO_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(GO_res_T2, gene_count.GO_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(GO_res_T3, gene_count.GO_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(GO_res_T4, gene_count.GO_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(GO_res_T5, gene_count.GO_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(GO_res_T6, gene_count.GO_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
target_ancestors <- c("GO:0001666","GO:0071456","GO:0002250","GO:0002376","GO:0002283","GO:0002252",
"GO:0002684","GO:0008152","GO:0030324","GO:0042990","GO:0006950",
"GO:0007610","GO:0008285","GO:0043522","GO:0010467","GO:0001501")
merged.res<-merged.res[merged.res$qvalue < 0.05,]
# Ancestor analysis
find_relation_to_ancestors = function (target_ancestors, GOs_to_check, ontology = 'BP')
{
if (nrow(GOs_to_check) == 0) {
return(data.frame(NULL))
}
GOs_to_check = GOs_to_check[GOs_to_check$ONTOLOGY == ontology,]
onto = switch(ontology, MF = GOMFANCESTOR, BP = GOBPANCESTOR,
CC = GOCCANCESTOR)
onto = AnnotationDbi::as.list(onto)
onto = stack(onto[GOs_to_check$ID])
relations_found <- data.frame(term_id = NULL, ancestor = NULL)
names_of_ancestor <- c()
for (GO in unique(onto$ind)) {
sub_onto <- onto[onto$ind == GO, ]
if (any(sub_onto$values %in% target_ancestors)) {
found_ancestor <- sub_onto$values[sub_onto$values %in%
target_ancestors][1]
temp_df <- data.frame(term_id = GO, ancestor = found_ancestor)
relations_found <- rbind(relations_found, temp_df)
}
else if (GO %in% target_ancestors) {
temp_df <- data.frame(term_id = GO, ancestor = GO)
relations_found <- rbind(relations_found, temp_df)
}
}
if (nrow(relations_found) == 0) {
message("NO terms associated with given ancestors were found")
return(data.frame(NULL))
}
num_ancestors <- length(unique(relations_found$ancestor))
cols <- RColorBrewer::brewer.pal(12, name = "Paired")
names(cols) <- unique(relations_found$ancestor)
relations_found$group_color <- cols[relations_found$ancestor]
GOs_ancestor_found <- GOs_to_check[GOs_to_check$ID %in%
relations_found$term_id, ]
#GOs_ancestor_found <- subset(GOs_ancestor_found, select = -c(group_color))
GOs_ancestor_found <- merge(GOs_ancestor_found, relations_found,
by.x="ID", by.y="term_id")
GOs_ancestor_found <- GOs_ancestor_found[order(GOs_ancestor_found$ancestor),
]
GOs_ancestor_found$group_color <- factor(GOs_ancestor_found$group_color,
levels = unique(GOs_ancestor_found$group_color))
goterms <- Term(target_ancestors)
goterms <- as.data.frame(goterms) %>% na.omit()
goterms$ancestor <- row.names(goterms)
colnames(goterms) <- c("ancestor_name", "ancestor")
GOs_ancestor_found <- merge(GOs_ancestor_found, goterms,
by = "ancestor")
return(GOs_ancestor_found)
}
GOs_ancestors_clust<-find_relation_to_ancestors(target_ancestors,merged.res)
num_cols<-length(unique(GOs_ancestors_clust$Cluster))
plt <-ggplot(GOs_ancestors_clust, aes(x = Cluster, y = Description, color=ancestor_name, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated'))
if(length(unique(GOs_ancestors_clust$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(color = guide_legend(override.aes = list(name = "", size = 8, shape = "\u25B2"))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
labs(shape="direction", colour="GO ancestor term") +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## GO gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
# Color by ancestor
plt
ggsave(filename = "gene_expression/figures/FigSn_GOenrichment_Abdominal.png", plt, width = 15, height = 23)KEGG abdominal
KEGG_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseKEGG_results.tab")
KEGG_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseKEGG_results.tab")
KEGG_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseKEGG_results.tab")
KEGG_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseKEGG_results.tab")
KEGG_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseKEGG_results.tab")
KEGG_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseKEGG_results.tab")
KEGG_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseKEGG_results.tab")
gene_count.KEGG_T0 <- KEGG_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T1 <- KEGG_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T2 <- KEGG_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T3 <- KEGG_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T4 <- KEGG_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T5 <- KEGG_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T6 <- KEGG_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(KEGG_res_T0, gene_count.KEGG_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(KEGG_res_T1, gene_count.KEGG_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(KEGG_res_T2, gene_count.KEGG_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(KEGG_res_T3, gene_count.KEGG_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(KEGG_res_T4, gene_count.KEGG_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(KEGG_res_T5, gene_count.KEGG_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(KEGG_res_T6, gene_count.KEGG_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
merged.res<-merged.res[merged.res$qvalue < 0.05,]
merged.res <- merged.res %>%
mutate(Description = factor(Description))
num_cols<-length(unique(merged.res$Cluster))
plt <-ggplot(merged.res, aes(x = Cluster, y = Description, color=type, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(name = "direction", values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated')) +
scale_color_manual(name = "direction",breaks = c('upregulated','downregulated'), values = c("#B31B21", "#1465AC"))
if(length(unique(merged.res$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## KEGG gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
plt
ggsave(filename = "gene_expression/figures/FigSn_KEGGenrichment_Abdominal.png", plt, width = 15, height = 23)MKEGG abdominal
MKEGG_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseMKEGG_results.tab")
MKEGG_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseMKEGG_results.tab")
MKEGG_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseMKEGG_results.tab")
MKEGG_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseMKEGG_results.tab")
MKEGG_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseMKEGG_results.tab")
#MKEGG_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseMKEGG_results.tab")
MKEGG_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseMKEGG_results.tab")
gene_count.MKEGG_T0 <- MKEGG_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T1 <- MKEGG_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T2 <- MKEGG_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T3 <- MKEGG_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T4 <- MKEGG_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
#gene_count.MKEGG_T5 <- MKEGG_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T6 <- MKEGG_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(MKEGG_res_T0, gene_count.MKEGG_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(MKEGG_res_T1, gene_count.MKEGG_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(MKEGG_res_T2, gene_count.MKEGG_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(MKEGG_res_T3, gene_count.MKEGG_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(MKEGG_res_T4, gene_count.MKEGG_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
#df_T5 <- left_join(MKEGG_res_T5, gene_count.MKEGG_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(MKEGG_res_T6, gene_count.MKEGG_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
#Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
merged.res<-merged.res[merged.res$qvalue < 0.05,]
merged.res <- merged.res %>%
mutate(Description = factor(Description))
num_cols<-length(unique(merged.res$Cluster))
plt <-ggplot(merged.res, aes(x = Cluster, y = Description, color=type, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(name = "direction", values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated')) +
scale_color_manual(name = "direction",breaks = c('upregulated','downregulated'), values = c("#B31B21", "#1465AC"))
if(length(unique(merged.res$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## MKEGG gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
plt
ggsave(filename = "gene_expression/figures/FigSn_MKEGGenrichment_Abdominal.png", plt, width = 13, height = 6)##WIKIPathway abdominal
WP_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseWP_results.tab")
WP_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseWP_results.tab")
WP_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseWP_results.tab")
WP_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseWP_results.tab")
WP_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseWP_results.tab")
WP_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseWP_results.tab")
WP_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Abdominal/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseWP_results.tab")
gene_count.WP_T0 <- WP_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T1 <- WP_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T2 <- WP_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T3 <- WP_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T4 <- WP_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T5 <- WP_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T6 <- WP_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(WP_res_T0, gene_count.WP_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(WP_res_T1, gene_count.WP_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(WP_res_T2, gene_count.WP_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(WP_res_T3, gene_count.WP_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(WP_res_T4, gene_count.WP_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(WP_res_T5, gene_count.WP_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(WP_res_T6, gene_count.WP_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
merged.res<-merged.res[merged.res$qvalue < 0.05,]
merged.res <- merged.res %>%
mutate(Description = factor(Description))
num_cols<-length(unique(merged.res$Cluster))
plt <-ggplot(merged.res, aes(x = Cluster, y = Description, color=type, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(name = "direction", values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated')) +
scale_color_manual(name = "direction",breaks = c('upregulated','downregulated'), values = c("#B31B21", "#1465AC"))
if(length(unique(merged.res$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## WP gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
plt
ggsave(filename = "gene_expression/figures/FigSn_WP_enrichment_Abdominal.png", plt, width = 13, height = 23)GO:BP PleuralPedal
GO_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseGO_results.tab")
GO_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseGO_results.tab")
GO_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseGO_results.tab")
GO_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseGO_results.tab")
GO_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseGO_results.tab")
GO_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseGO_results.tab")
GO_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseGO_results.tab")
gene_count.GO_T0 <- GO_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T1 <- GO_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T2 <- GO_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T3 <- GO_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T4 <- GO_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T5 <- GO_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.GO_T6 <- GO_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(GO_res_T0, gene_count.GO_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(GO_res_T1, gene_count.GO_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(GO_res_T2, gene_count.GO_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(GO_res_T3, gene_count.GO_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(GO_res_T4, gene_count.GO_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(GO_res_T5, gene_count.GO_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(GO_res_T6, gene_count.GO_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
target_ancestors <- c("GO:0001666","GO:0071456","GO:0002250","GO:0002376","GO:0002283","GO:0002252",
"GO:0002684","GO:0008152","GO:0030324","GO:0042990","GO:0006950",
"GO:0007610","GO:0008285","GO:0043522","GO:0010467","GO:0001501")
merged.res<-merged.res[merged.res$qvalue < 0.05,]
GOs_ancestors_clust<-find_relation_to_ancestors(target_ancestors,merged.res)
num_cols<-length(unique(GOs_ancestors_clust$Cluster))
plt <-ggplot(GOs_ancestors_clust, aes(x = Cluster, y = Description, color=ancestor_name, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated'))
if(length(unique(GOs_ancestors_clust$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(color = guide_legend(override.aes = list(name = "", size = 8, shape = "\u25B2"))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
labs(shape="direction", colour="GO ancestor term") +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## GO gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
# Color by ancestor
plt
ggsave(filename = "gene_expression/figures/FigSn_GOenrichment_Pleural.png", plt, width = 15, height = 30)KEGG_pleural/pedal
KEGG_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseKEGG_results.tab")
KEGG_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseKEGG_results.tab")
KEGG_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseKEGG_results.tab")
KEGG_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseKEGG_results.tab")
KEGG_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseKEGG_results.tab")
KEGG_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseKEGG_results.tab")
KEGG_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseKEGG_results.tab")
gene_count.KEGG_T0 <- KEGG_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T1 <- KEGG_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T2 <- KEGG_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T3 <- KEGG_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T4 <- KEGG_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T5 <- KEGG_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.KEGG_T6 <- KEGG_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(KEGG_res_T0, gene_count.KEGG_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(KEGG_res_T1, gene_count.KEGG_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(KEGG_res_T2, gene_count.KEGG_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(KEGG_res_T3, gene_count.KEGG_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(KEGG_res_T4, gene_count.KEGG_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(KEGG_res_T5, gene_count.KEGG_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(KEGG_res_T6, gene_count.KEGG_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
merged.res<-merged.res[merged.res$qvalue < 0.05,]
merged.res <- merged.res %>%
mutate(Description = factor(Description))
num_cols<-length(unique(merged.res$Cluster))
plt <-ggplot(merged.res, aes(x = Cluster, y = Description, color=type, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(name = "direction", values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated')) +
scale_color_manual(name = "direction",breaks = c('upregulated','downregulated'), values = c("#B31B21", "#1465AC"))
if(length(unique(merged.res$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## KEGG gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
plt
ggsave(filename = "gene_expression/figures/FigSn_KEGGenrichment_Pleural.png", plt, width = 15, height = 20)MKEGG pleural/pedal
MKEGG_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseMKEGG_results.tab")
MKEGG_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseMKEGG_results.tab")
MKEGG_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseMKEGG_results.tab")
MKEGG_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseMKEGG_results.tab")
MKEGG_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseMKEGG_results.tab")
MKEGG_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseMKEGG_results.tab")
MKEGG_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseMKEGG_results.tab")
gene_count.MKEGG_T0 <- MKEGG_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T1 <- MKEGG_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T2 <- MKEGG_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T3 <- MKEGG_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T4 <- MKEGG_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T5 <- MKEGG_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.MKEGG_T6 <- MKEGG_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(MKEGG_res_T0, gene_count.MKEGG_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(MKEGG_res_T1, gene_count.MKEGG_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(MKEGG_res_T2, gene_count.MKEGG_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(MKEGG_res_T3, gene_count.MKEGG_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(MKEGG_res_T4, gene_count.MKEGG_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(MKEGG_res_T5, gene_count.MKEGG_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(MKEGG_res_T6, gene_count.MKEGG_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
merged.res<-merged.res[merged.res$qvalue < 0.05,]
merged.res <- merged.res %>%
mutate(Description = factor(Description))
num_cols<-length(unique(merged.res$Cluster))
plt <-ggplot(merged.res, aes(x = Cluster, y = Description, color=type, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(name = "direction", values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated')) +
scale_color_manual(name = "direction",breaks = c('upregulated','downregulated'), values = c("#B31B21", "#1465AC"))
if(length(unique(merged.res$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## MKEGG gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
plt
ggsave(filename = "gene_expression/figures/FigSn_MKEGGenrichment_pleural.png", plt, width = 15, height = 6)##WIKIPathway pleural
WP_res_T0 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_0/ gseWP_results.tab")
WP_res_T1 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_1/ gseWP_results.tab")
WP_res_T2 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_2/ gseWP_results.tab")
WP_res_T3 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_3/ gseWP_results.tab")
WP_res_T4 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_4/ gseWP_results.tab")
WP_res_T5 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_5/ gseWP_results.tab")
WP_res_T6 <- read_delim(file = "gene_expression/output/TS_output/LCvsWC_Pleural/GSEA_results/Lab_cross_vs_Wild_cross_TP_6/ gseWP_results.tab")
gene_count.WP_T0 <- WP_res_T0 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T1 <- WP_res_T1 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T2 <- WP_res_T2 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T3 <- WP_res_T3 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T4 <- WP_res_T4 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T5 <- WP_res_T5 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
gene_count.WP_T6 <- WP_res_T6 %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
df_T0 <- left_join(WP_res_T0, gene_count.WP_T0, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T1 <- left_join(WP_res_T1, gene_count.WP_T1, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T2 <- left_join(WP_res_T2, gene_count.WP_T2, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T3 <- left_join(WP_res_T3, gene_count.WP_T3, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T4 <- left_join(WP_res_T4, gene_count.WP_T4, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T5 <- left_join(WP_res_T5, gene_count.WP_T5, by = "ID") %>% mutate(GeneRatio = count/setSize)
df_T6 <- left_join(WP_res_T6, gene_count.WP_T6, by = "ID") %>% mutate(GeneRatio = count/setSize)
merged.res <- as.data.frame(merge_result(list(Ctrol=df_T0,
Hyp_2h=df_T1,
Hyp_6h=df_T2,
Hyp_5d=df_T4,
Hyp_6d=df_T5,
Reox=df_T3,
Reco=df_T6)))
## Set up/downregulation
merged.res$type = "upregulated"
merged.res$type[merged.res$NES < 0] = "downregulated"
merged.res<-merged.res[merged.res$qvalue < 0.05,]
merged.res <- merged.res %>%
mutate(Description = factor(Description))
num_cols<-length(unique(merged.res$Cluster))
plt <-ggplot(merged.res, aes(x = Cluster, y = Description, color=type, shape=type)) +
geom_point(aes(size=setSize,shape=type),show.legend = TRUE) +
scale_shape_manual(name = "direction", values=c("\u25B2","\u25BC"),breaks = c('upregulated','downregulated')) +
scale_color_manual(name = "direction",breaks = c('upregulated','downregulated'), values = c("#B31B21", "#1465AC"))
if(length(unique(merged.res$setSize))>1){
plt<-plt+ scale_size(name = "term size",
range = c(5,15))
}
plt<-plt+
scale_y_discrete(labels = function(x) str_wrap(str_replace_all(x, "foo" , "\n"),width = 60)) +
theme_bw(base_size=15) +
theme(axis.text.x=element_text(angle=-45,vjust=0)) +
guides(shape = guide_legend(override.aes = list(size = 8))) +
guides(size = guide_legend(override.aes = list(shape = "\u25B2"))) +
ylab("") +
xlab("") +
facet_wrap(~Cluster,scales='free_x',ncol=num_cols)
## WP gene set enrichment
# Upregulated (Higher expression in Lab_cross)
# Downregulated (Lower expression in Lab_cross)
plt
ggsave(filename = "gene_expression/figures/FigSn_WP_enrichment_Pleural.png", plt, width = 13, height = 23)